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ABSTRACT 

The massive stars in the Galactic center inner arcsecond share analogous properties with the so-called Hot 
Jupiters. Most of these young stars have highly eccentric orbits, and were probably not formed in-situ. It has 
been proposed that these stars acquired their current orbits from the tidal disruption of compact massive binaries 
scattered toward the proximity of the central supermassive black hole. Assuming a binary star formed in a thin 
gaseous disk beyond 0.1 pc from the central object, we investigate the relevance of disk-satellite interactions 
to harden the binding energy of the binary, and to drive its inward migration. A massive, equal-mass binary 
star is found to become more tightly wound as it migrates inwards toward the central black hole. The migration 
timescale is very similar to that of a single-star satellite of the same mass. The binary's hardening is caused 
by the formation of spiral tails lagging the stars inside the binary's Hill radius. We show that the hardening 
timescale is mostly determined by the mass of gas inside the binary's Hill radius, and that it is much shorter 
than the migration timescale. We discuss some implications of the binary's hardening process. When the 
more massive (primary) components of close binaries eject most their mass through supernova explosion, their 
secondary stars may attain a range of eccentricities and inclinations. Such processes may provide an alternative 
unified scenario for the origin of the kinematic properties of the central cluster and S-stars in the Galactic center 
as well as the high velocity stars in the Galactic halo. 

Subject headings: accretion, accretion disks — binaries: general — Galaxy: center — hydrodynamics — 
methods: numerical 



1. INTRODUCTION 

The supermassive black hole (SMBH) in the Galactic center 
is surrounded by a parsec-scale cluster of young and massive 
stars. At separations larger than ~ 1" or 0.04 pc , about 100 
OB stars are observed in a moderately thin disk that rotates 
clockwise on the sky and extends up to ~ 0.6 pc. The disk's 
scale height to radius ratio is ^0.1. The stars have masses 
> 10 Mg and a common age ~ 6 Myr. Thei r eccentricities 
range up to 0.8, with a typ ical value ~ 0.4 (e.g jLu et al.l2009l: 
iBartko etai . 2009, 2010). There is observational evidence for 
the existence of a second, less massive counter-clockwise disk 
comprising stars with si milar age and kinema t ic properties 
as in the cl ockwise disk dPaumard et al.l 12006: Bart ko et al.l 
120091 120101) . In contrast, stars within the central arcsecond 
have orbits with random orientations and even higher eccen- 
tricities. Known as the "S-stars", they are main-sequence B 
stars with masses of ~ 10 and a lifespan of several x 10 
yrs (e .g., ISchodel et ai]|2002t iGhez et al.ll200l IBartko et al l 
2010). The orbits of abo ut 20 of them have been determined 
bv lGillessen et al.l J2009). Their semi-major axes are as small 
as ~ 4 x 10~ 3 pc 800 AU, and their eccentricities are typi- 
cally > 0.8. 

Explaining the presence of such young stars so close to 
the SMBH is challenging. The standard star formation sce- 
nario, invoking the collapse of self-gravitating cold molecular 
clouds, is unlikely to occur , because of t he strong tidal per- 
turbation from the SMBH (lMorrislll993l) . Nonetheless, the 
tidal disruption of a molecular cloud can form a thin a ccretion 
disk around the central object (Bonnell & Rice 2008). In-situ 



star formation may then be triggered when the disk density is 
high enoug h so its self-gravity overcomes the tidal force of the 
SMBH dLevin & Beloborodovl 120031: iNavakshin & Cuadral 
120051) . a process first predi cted for the massive accretion 
disks present in AGNs (e.g., | Kolvkhalov & Svunvaev|[r980l : 
IShlosman & Begelmadll989HGoodmanll2003l) 

Alternatively, stars may have formed far from their current 
location and have migrated inwards due to dynamical friction. 
While the migration time for individual stars is much too long, 
a massive (> 10 5 M Q ), dense stellar cluster formed several 
parse cs away could reach the central parsec in a few 10 6 yrs 
(e.g.. lGerha rd 2001). Both the self-gravitating disk scenario 
and the infalling c luster scenario can produce stellar disks 
(see lYu et alj|2007l and references therein), but the observed 
properties of the stellar disk in the central parsec favor in-situ 
formatio n (Navakshi n & Sunv aev 2005: lPaumard et aT]|2006l : 
ILu et alj|2009l: IBartko et alj|2009l) . 

Additional complexity arises when attempting to reproduce 
the peculiar orbits of the S-stars, which are even closer to the 
SMBH. The most promising mechanism involve s the three - 
body interacti on of a binary star and th e SMBH dHillsll 1 9 8 8h . 
In particular, iGould & Quilled d2003l) proposed that the S- 
stars were initially in compact binary systems. The binary 
can get tidally disrupted if it is scattered sufficiently close to 
the SMBH, such that its binary separation is a significant frac- 
tion of its Hill radius at its closest approach to the SMBH. In 
addition to potentially explaining the orbits of the S-stars, the 
binary's tidal disruption s cenario cou l d account for the ob- 
served hypervelocity s t ars dHil ls 1988; Yu & Tremaindl2003l : 
iGualandris e"tal]l2005l: IL~u et al.H2007tlPeretsll2009l). most of 
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iicojiWMW 1 are B-type stars, just like the S-stars (e.g., Brown et al. 
20071) . In this context. lLu et al.l (2010) and Zhang et al.ld201( 
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have recently pointed out that the spatial distribution of the 
hypervelocity stars is consistent with being located on two 
thin disks. These authors have shown that this distribution 
could be reproduced by assuming that the hypervelocity stars 
originate from the tidal disruption of tight binaries initially 
located on two stellar disks near the Galactic center. In the 
tidal disruption scenario, one component of a compact bi- 
nary is ejected at a velocity comparable to its orbital velocity 
around the binary's center of mass, while the other compo- 
nent remains bound to the central black hole, with eccentric- 
ity and inclination very close to those of the binary's center 
of mass before disruption. To produce bound B-type stars 
with eccentricities comparable to those of the S-stars, the 
original separation of the incoming binary (d istance between 
both stars) would have to be ~ 1 AU (e.g.. lGould & O uillen 
2003). Producing unbound B-type stars with ejection veloci- 
ties ~ 2000 km s would demand an original sep aration typ- 
ically < 0. 1 AU (iHillsll 19881: iBromley et al]l2006l) . 

The presence of a couple of compact b inary stars has been 
confirmed in the Galactic center (e.g.. [Martins et all 120061) . 
Several processes may account for their formation. The dy- 
namic relaxation of a population of embedded single stars 
may lead to their capture into binary-star systems. Alterna- 
tively, binary stars may form during the fra gmentation of a 
self-gravitating disk (Alexand er et al.l l2008a). But, these new 
born binaries generally have wide separations, and they are 
easily disrupted at large distance from the central black hole. 
The ejection speed of such disrupted binaries would generally 
be too small to produce either the hypervelocity stars or the 
large eccentricities of the S-stars. To do so, additional pro- 
cesses that can harden the binaries and drive their inward mi- 
gration should be at work. Several stellar dynamical processes 
have been considered to bring the necessary binaries to the 
proximity of the SMBH (|Perets et alj|2007t Lockman n et al.l 
120081; iMadigan et al.ll2009t iPerets & Gualandrisll2010l) . 

Interestingly, the S-stars have many properties analogous 
to the observed exoplanets, in particular the so-called Hot 
Jupiters. The latter are planets typically as massive as Jupiter, 
orbiting within 0.5 AU from their host star. S ome of these 
planets have eccentricities up to 0.97 (see e.g.. ISantosll2008L 
for a recent review on exoplanets properties). According to 
planet formation theories, it is very unlikely to for m a giant 
gaseo us planet at such small distances from a star dLin et al.l 
1996). Instead, giant planets are believed to form at larger 
separations in the protoplanetary gaseous disk surrounding 
the central object. The tidal interaction between the disk and 
the p lanet genera lly leads to decrease the planet's semi-major 
axis (IWardlll997h . This is known as planetary migration. The 
migration of low-mass planets (say, an Earth-mass planet or- 
biting a solar mass star) has received considerable attention 
as it may occur on timescales much shorter than the life- 
time of the protoplanetary disk (Tana kaet al.l 120021) . poten- 
tially bringing all these planets to the close proximity of their 
central host. 

In the central parsec of our Galaxy, the star-to-supermassive 
black hole mass ratio is close to the Earth-to-Sun mass ra- 
tio. It suggests that "planet-like migration" could be an 
effici ent way to bring young stars to the inner arc secon d 
(g^Artvmow icz et al.l 11991 iLin et al.|[l99H IChangl 2008). 
iLevin " (120071) has proposed a scenario in which the S-stars 
formed in a sub-pc self-gravitating disk and migrated in- 



wards by interacting with the disk 2 . Because disk-satellite 
interactions damp the satellite's eccentricities and inclina- 
tions on short timescales, an additional process is required 
to explain the large dispersion of the S- stars eccent ricities 
and inclinations, like for the exoplanets. iLevinl (120071) pro- 
posed a fast resonant relaxation process th at would randomize 
the S -stars' eccentricity and inclination dRauch & Tremainel 
[19961) . The presence of an intermediate-mass black hole may 
also influence stellar orbi ts, increasing their eccentricity and 
inclination considerably (|Yu et all l2007h iMerritt et al.l 120091 
iGualandris & Merrittll2009l) . 

Our approach is to investigate the tidal interaction between 
a massive binary star and its parent gaseous disk as a mecha- 
nism to harden the binary, and to bring it to the SMBH prox- 
imity. This migration scenario differs from the stellar scatter- 
ing scenario, where stars have parabolic encounters with the 
central black hole. The outline of this paper is the following. 
The physical properties of the system we model are described 
in § 2] We then discuss in § [3| the migration timescales we 
expect for a single star. The numerical set-up used in our hy- 
drodynamical simulations is detailed in § [4] To alleviate the 
resolution requirements, we show in §[5] that it is possible to 
treat an equivalent rescaled problem. The calculation results 
obtained with the rescaled problem are presented in § We 
then come back to our original problem in § [7] and we finally 
discuss our results in the context of the massive stars near the 
Galactic center in § [8] 

2. PHYSICAL MODEL 

The system that we study comprises a supermassive black 
hole of mass M sm bh = 3 x 10 6 M Q surrounded by a gaseous 
disk, where stars are assumed to have already formed by disk 
fragmentation. For simplicity, only one satellite, a binary star, 
is embedded in the disk. Our aim is to investigate how disk- 
satellite interactions alter the orbital elements of the binary 
star (orbital migration, and possible evolution of the distance 
between both stars). The case of a single-star satellite with the 
same mass is also considered for comparison. The p roper ties 
o f the binary star and of the disk are described in § 12. 1 1 and 

§1121 

2.1. Binary properties 

We consider a binary with equal-mass stars. Both stars 
have a fixed mass = 15 Mp>, akin to the star SO-2 or- 
biting Sgr A* (Mart ins et alj|2008l) . Gas accretion onto the 
stars is discarded. The satellite to primary mass ratio is 
q = 2M*/M sm bh ~ 10~ 5 , which is about three times the Earth- 
to-Sun mass ratio. The motion of the stars with respect to 
the binary's center of mass is prograde. Resu lts with a retro- 
grade binary will also be presented in § 16.41 We denote by 
a and e the semi-major axis and eccentricity of the binary's 
center of mass, with respect to the central black hole. The 
semi-major axis and eccentricity of the binary, with respect 
to its center of mass, are denoted as abin an d e bin- We assume 
that the binary has initially e = 0, a = 0.1 pc, ebin = 0, and 
flbin = 5 x 10~ 4 pc s» 100 AU. The binary's initial semi-major 
axis is thus «bin = 0.3Rn, where Rn = a is the bi- 

nary's Hill radius. 

2.2. Properties of the gaseous disk 



2 Not necessarily the same disk that originated the OB stars observed now. 
The S-stars could have formed in a previous star formation event. 
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The properties of gaseous disks around massive black holes 
prior to star formation have been descr ibed in many studies 
(e.g., En & Pringle][l987t lLevinll2007h . Such disks are ex- 
pected to be geometrically thin, with a typical aspect ratio 
h (ratio of the disk's scale height H to orbital ra dius r) be- 
tween 0.5% and 1% (iNayakshin & Cuadrall2005h . But, in- 
vestigating star-disk interactions requires knowledge of the 
disk properties after star formation, which are largely uncer- 
tain. The decrease of the disk density, and the increase of its 
temperature depend on the number and on the mass of stars 
formed (e.g., [Sh losman & Begelman 1989; Nayakshinl l2006c 
Nava kshin et al.l 120071) . We assume that the increase of the 
disk temperature after star formation is mild, and we set the 
disk aspect ratio to h = 1% at r = 0.1 pc (the initial loca- 
tion of our binary-star system). This yields Rh w 1.5H at the 
binary's initial location, which makes the two-dimensional 
approximat ion for the disk ap plicable in absence of gas ac- 
cretion (D'Angelo et al. 2005). In addition, the fact R H > H 
indicates that, despite its rather small mass, the binary may 
open a gap around its orbit, the depth of which depends on the 
disk's turbulent viscosity around the binary's orbit. The lat- 
ter is modelled by a kinematic viscosity v. The gap-opening 
criter ion takes the form dLin & PapaloizouH l993; Cri da et all 
2006) 

3 Mo) m-" + Jo^L <,. 

4 \3/ qa L Ll(a) 

For our disk and binary parameters, the first term in Eq. (fl} 
is s» 0.5. Assuming v = ah 2 a 2 £l (IShakura & Sunvaevl 1973b . 
the second term in Eq. (fl} equals 50a(a)h 2 (a)/q. Values of a 
in gravito-turbulent, critically fragmenting disks should typi- 
cally range from 10~ 2 to IP" 1 (e.g.. iLin et alJI 1988c [Gammie 
I200U iLodato&RicelllOOlt lRafikovll2009lh Still, the value 
of a after star formation is uncertain, as stars modify the 
disk's local properties (heating by shocks, gap-opening). For 
a (a) = 10~ 2 , the second term in Eq. (HJ equals 5, suggesting 
the binary would slowly open a shallow gap. But, note that 
its actual value depends sensitively on h(a). Given the un- 
certainties on both h and a, we assume that the viscosity at 
the binary's location is low enough so that the binary rapidly 
clears a clean gap around its orbit. For this purpose, we take 
a(a) = 10~ 3 , and the second term in Eq. (Q]i equals 0.5. 

We now describe the properties of the unperturbed disk 
before the satellite opens a gap. The disk's angular veloc- 
ity D. is slightly sub-Keplerian, the gravitational force due to 
the central object being compensated for by the centrifugal 
force and by the pressure force. The temperature T and sur- 
face density E are set as power-law functions of r: T(r) — 
7().ipc('"/0.1pc) _1 , andE(r) = E .i pc (r/0.1 pc) -1 / 2 . Assum- 
ing a mean molecular weight jj. — 2.4 and an adiabatic index 
7=5/3, the initial temperature at O.lpc is 7o.i pc ~2x 10 3 K. 
The initial surface density at 0.1 pc, Eo.i pc . is a free parameter 
in our study. It ranges from ~ 20 g cm~ 2 to ~ 200 g cm -2 for 
Go.ipc between 1 and 10, where Qo.ipc stands for the Toomre 
parameter at 0.1 pc, 



with c s the local sound speed. For simplicity, the gas self- 
gravity is neglected in our study. 

If stars formed along the lines of the gravitational instabil- 
ity's model, the disk is expected to be radiatively efficient in 



the stars forming region. After star formation, the increase of 
the gas temperature due to stellar heating, and the decrease 
of the gas density c ould significantly decrease the disk's opti- 
cal thickness (e.g.. lBell & Linl[l994l) . It is therefore possible 
that the gas becomes radiatively inefficient after star forma- 
tion. Given again the uncertainties on the disk properties af- 
ter star formation, we assume for simplicity that the gaseous 
disk is radiatively efficient, modeled with a locally isothermal 
equation of state (the initial temperature profile is maintained 
constant with time). 

3. MIGRATION TIMESCALE OF A SINGLE STAR 

We briefly review in this section the different regimes of mi- 
gration resulting from the tidal interaction between a gaseous 
disk and a single-star satellite. A d etailed review on disk- 
satellite interactions can be found in Masset (120081) . Migra- 
tion timescales are given assuming a two-dimensional disk 
with surface density and temperature profiles decreasing as 
r~ a and r~P, respectively. The torque T exerted by the disk 
on the satellite alters the satellite's semi-major axis and ec- 
centricity. Assuming the satellite is initially on a circular or- 
bit, the time evolution of its semi-major axis a is given by 
2B<7M sm bh a a = F, where B is the second Oort's constant. 
The satellite's migration timescale T is defined as T = \a/d\. 

3.1. Type I migration 

Type I migration applies to satellites that do not clear a 
gap around their orbit. The gap-opening criterion is recalled 
in Eq. £[). Assuming a di sk aspect ratio at the star's lo- 
cation of ~ 1% (see § I2.2K and a central object of a few 
10 6 solar masses, type I migration typically applies to stars 
up to a few solar masses. In this migration regime, the 
torque Tj exerted by the disk on the satellite is the sum of 
the differential Lindblad torque a nd of the horseshoe drag 
(Paardekooper & Papaloizou 2009). The former accounts for 
the angular momentum carried away by the spiral density 
waves generated by the satellite, whereas the latter corre- 
sponds to the angular momentum exchanged between the 
satellite and the gas inside its corotation region (also known 
as horseshoe region). The total torque can be written as 

Tj = -cVEa 4 £l 2 /r 2 , (3) 

where all the unperturbed disk's quantities in Eq. (f3j) are 
to be evaluated at the satellite location. The dimension- 
less factor Q depends on a, j3, on the softening length 
of the satellite's gravitational potential (required to adjust 
semi-analytic estimates of the horseshoe drag with the re- 
sults of two-dimensional hydrodynamical simulations). It 
also depends on the properties of the gas dis k, among which 
its self-gravity (Baruteau & Masset 2008b), turbulent vis- 
cosity (Masset 2001; Barute au & Linl 120101). and radiative 
properties dBaruteau & Massetl2008al;lMasset & Casolil2009c 
Paardekooper et al. 2010a b). In a two-di mens ional, non self- 
gravitating locally isothermal disk (see § 12.2b . C\ can be ap- 
proximated 3 as C\ « 0.8 + 1.0<7 + 0.9/3 dPaardekooper et all 
l2010ai their equation 49). For any density and temperature 
profiles decreasing with radius, Q is positive, and type I mi- 
gration is directed inwards. Using previous notations, the 

3 The softening length to scale height ratio is taken equal to 0.4, and the 
horseshoe drag is assumed to be fully unsaturated. For more details on the 
horseshoe drag evaluat ion in isothermal disk models, the reader is referred 
e.g. to lBaruteau & Linl 120101) . 
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timescale for type I migration, Tj, reads 
Zi « 4.3 x 10 7 yrsx (0.8 + 1.0<r + 0.9/3) 



-l / « 



3 x 10- 6 



10 



-2 



Go. 



1 pc 



30 



M, 



smbh 



3 x 10 6 M F 



-l 



(4) 



,0.1 pc 

In the particular case where ct = 2 and j8 = 1, Eq. has 
the same dependence with the disk and satellite parameters 
as in equation (40) of lLevinl (120071) . With our model param- 
eters (ct = 1/2, P = 1, h = 10" 2 , M smbh = 3 x 10 6 M Q and 
a = 0.1 pc), and further assuming q = 10~ 7 (/?h ps 0.3//) and 
Go.i pc = 3, we find Z\ ~ 6 x 10 7 yrs ps 3.5 x 10 4 orbital peri- 
ods at 0.1 pc. This timescale is slightly longer than the life- 
time of the stars near the Galactic center. Note that ct = 1/2 
and j3 = 1 make Ti independent of a. With ct = 2 and j8 = 1, 
Tj oc a 3 / 2 and goes down to a few 10 6 yrs at a ~ 0.01 pc. Note 
also that, since Tj « q~ l , more massive stars (but still sub- 
ject to type I migration) will migrate faster. However, based 
on our thin disk assumption, the S-stars would open a gap 
around their orbit. They are therefore more likely subject to 
the slower type II migration. Still, it is possible that the S- 
stars have substantially migrated in before they acquired their 
asymptotic mass. 

3.2. Type II migration 

Type II migration concerns satellites that are massive 
enough to progressively deplete their horseshoe region, 
thereby opening a gap around their orbit. Assuming again 
a disk aspect ratio of ~ 1%, and a central object of a few 
10 6 Mq, type II migration typically applies to stars more mas- 
sive than ~ 10 M©. Compared to the type I migration regime, 
the horseshoe drag is much reduced, and the Lindblad torque 
balances the viscous torque exerted by the disk. The to- 
tal torque on the satellite can be written as a fraction Cn of 
the vi scous torque due to the outer disk (Crida & Morbidelli 
2007, their equation 15). The factor C\\ features the time- 
dependent fraction of gas / gas left in the satellite's horseshoe 
region. 

The particular case where / gas tends to zero is usually called 
the standard type II migration regime. Its timescale reads 



Hi ' 



2ri 



3v(r ) 



47T£(r )r 2 1 



(5) 



where M sat = gM sm t,h denotes the satellite mass, and where E 
is the density of the disk perturbed by the satellite. In Eq. (0, 
r a is the location in the outer disk where most of the satel- 
lite's angular momentum is deposited. It can be approximated 
as the location of the outer separatrix of the satellite's horse- 
shoe region, given by r p ps a + 2.5Rn for gap-opening satel- 
lites d^iMetil]|2Sl). 

When the satellite mass is much smaller than the local disk 
mass 47r£(r )r 2 , the type II migrati on timescale correspond s 
to the disk's local viscous timescale dLin & Pa paloizou 1986). 
In this regime, known as disk-dominated type II migration, the 
migration timescale Tn d can be written as 



Tn,d ~ 1 • 



xloVsx^) 



-1 



'smbh 



3 x 10 6 M. 



lO- 2 

3/2 



-2 



0.1 pc 
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where a and h are to be evaluated at r D . This timescale is 
typically one order of magnitude longer than the timescale 
for type I migration, given in Eq. 

In the opposite case where the satellite mass is large com- 
pared to the local disk mass (satellite-dominated type II mi- 
gration), the migration timescale Tn )S is 



Tli, s « Tn d x 



/ M Si 



V47r£(r )r 2 
where the satellite to local disk mass ratio is 



(7) 



47r£(r )r 2 



1.7 x 10 



\l0Me 



10 2 g cm 



-l 



0.1 pc 



-2 



(8) 



Eq. © shows that in the early stages of their formation and 
evolution, most of the massive stars near the Galactic center 
should be subject to the disk-dominated type II migration, and 
thus migrate in about the disk's viscous timescale. Depletion 
of the gas disk, or substantial migration towards the central 
object, should however slow down migration as the stars iner- 
tial mass becomes comparable to the local disk mass. 

3.3. Type III migration 

Migrating satellites that open a partial gap experience an 
additional corotation torque due to fluid elements flowing 
across the horseshoe region (moving from the inner disk to 
the outer disk if the satellite migrates inwards). At small 
migration rates a, this additional corotatio n torque is pro- 
portional to a (Masset & Papaloizou 2003). When its am- 
plitude is large enough, this torque may thus trigger a run- 
away migration. Runaway occurs, roughly speaking, when 
the mass deficit of the horseshoe region (the difference be- 
tween the mass of the horseshoe region, and the mass that it 
would have if it had a uniform surface density equal to the 
density of the orbit-crossing flow) exceeds the satellite's mass 
dMasset & P apaloizou 2003). The occurrence and timescale 
of runaway migration, also referred to as type III migra- 
tion, depend sensitively on the parameters entering the gap- 
opening criterion ( q, h and a at the satellite's location), and 
on the disk mass dMasset & Papaloizou! 120031) . It concerns 
intermediate-mass satellites (massive satellites clear a wide, 
deep gap, and the density of the orbit-crossing flow is too 
small to induce a significant mass deficit) in massive disks 
(the more massive the disk, the larger the density of the 
orbit-crossing flow). This migration regime could then be 
particularly relevant to > 3 — 10 Mr., stars embedded in thin 
(h ~ 1%) massive (Q < 10) disks. Numerical simulations of 
disk-satellite interactions show that, depending on the resolu- 
tion of the gas flow surrounding the satellite, the timescale for 
runaway migration may be as short as a few 10 2 orbits (e.g., 
iMasset & Papaloizoull2003l: ICrida et alj|2009l) . It is thus pos- 
sible that some of the massive stars near the Galactic center 
formed at separations > 0.1 pc, and reached the proximity of 
the central black hole through type III runaway migration. 

4. NUMERICAL MODEL 

4.1. Numerical method 

All our calculations are performed with the two- 
dimensional code FARGO. It is a staggered mesh code that 
solves the Navier-Stokes and continuity equations on a polar 
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grid. An upwind transport scheme is used along with a har- 
monic, second-order slope limiter (Ivan LeeJ[T977h . Its main 
feature is to use a change of rotating frame on each ring of 
the g rid, which increases the timestep significantly dMassell 
12000) . 

Our calculation results are expressed in the following units. 
The initial semi-major axis oq of the binary's center of mass 
(or of the satellite if single) is the length unit, the mass of the 

central objectM sm b n is the mass unit, and (GM sm bh/flo)~ 1,/2 is 
the time unit, G being the gravitational constant (G = 1 in our 
unit system). Whenever time is expressed in orbital periods, 
it refers to as the orbital period of the binary's center of mass 
(or of the satellite if single) around the central black hole. 

We describe in this paragraph the numerical set-up of our 
calculations. Unless otherwise stated, the disk and satellite 
parameters are those described in § [2] Since our model as- 
sumes a locally isothermal disk, simulations do not include an 
energy equation (the temperature's initial radial profile does 
not evolve with time). All calculations are performed in the 
frame corotating with the binary's center of mass, or with the 
satellite if single. Wave-killing zones are used next to the 
grid's inner and ou ter edges in order to min imize unphysi- 
cal wave reflexions (Ide Va l-Borro et al. 2006 ). In simulations 
with a binary satellite, we follow ICresswell & Nelsonl ([2006) 
and set the timestep as the minimum value between the hy- 
drodynamical timestep, and a fixed fraction of the binary's 
internal orbital period, 



St] 



bin 



In 
400 



4n 



2GM* 



1/2 



(9) 



where d^ m stands for the separation between the two stars. To 
avoid a violent initial evolution of the disk, the stars' gravi- 
tational potential is slowly turned on over 10 orbital periods. 
Discussion on the for m of the stars potential, and on the grid 
resolution follows in 8 14.21 



4.2. Resolution issues 

The time-evolution of the binary should sensitively depend 
on how the gas flow around and between the stars is resolved. 
As is usual in disk-satellite 2D hydrodynamical simulations, 
the stars potential is smoothed over a softening length £. Two 
parameters therefore determine the resolution of our calcu- 
lations: the grid's resolution and the softening length of the 
stars potential. In this study, we choose a rather large value 
for e in order to prevent a large accumulation of gas around 
each star, which would otherwise require a prohibitively high 
grid's resolution to be properly handled. The default value for 
e is half the binary's initial semi-major axis abin.o (the same 
value is used in simulations with a single satellite). The soft- 
ening length is therefore kept constant througho ut th e simula- 
tion. Smaller values of e will be considered in § 16.31 

We aim to resolve the satellite's Hill radius with about 50 
cells along the radial and azimuthal directions (ICrida et al.l 
2009). This requirement implies a maximum size for the 
grid cells equal to 5 max = 0.02R H ~ 3 x 10~ 4 a. To model the 
global disk-satellite interaction, it is enough to set the grid's 
radial extent to Ar = ISR^. This sets the number of cells along 
the radial direction to N r = 750. The cell number along the az- 
imuthal direction is then N s = 2 7£a/ 5 max « 21 000. In contrast 
to N T , the value of N s is about one order of magnitude larger 
than its typical value in 2D hydrodynamical simulations of 
disk-satellite interactions. To maintain a reasonable compu- 
tational cost, we show in § [5] that it is possible to transform 



Table 1 

Parameters of the calculations of §[5] where the rescaling method is 

illustrated. All runs have the same values for d[ = q/h 3 and 
do = qa 2 Q./v, the dimensionless quantities entering the gap-opening 
criterion. 



Parameter 


Run a 


Run b 


Run c 


h 


1% 


3% 


5% 


q<xh? 


8 x 1(T 6 


2.2 x 10~ 4 


io- 3 


v/a 2 £l « q 


8 x 1(T 8 


2.2 x 10~ 6 


10- 5 


a = v/h 2 a 2 Q. 


8 x 1(T 4 


2.4 x 10~ 3 


4x 10~ 3 


Ryi/a h 


1.4% 


4.2% 


6.9% 


^mm i '"max 


0.88,1.16 


0.64,1.48 


0.4,1.8 


N T ,N S 


448,6700 


448,2238 


448,1344 



our problem into an equivalent rescaled problem, less com- 
putationally demanding, by making use of the dimensionless 
parameters in the gap-opening criterion at Eq. (Q]). 

5. RESCALED PROBLEM 

For a single satellite with fixed semi-major axis a, the 
gap-opening criterion involves two dimensionless quantities: 
d\ = q/h? and di = qa 2 £l/v. Any changes of q, h and v 
that keep d\ and di constant will preserve the gap's density 
profile (width, depth) in scaled units A/Ru, where A stands 
for the disk-satellite radial distance ( C rida et al. 2006j). The 
major resolution issue pointed out in § 14.21 is due to the fact 
N s oc R H _1 oc q ~ ! / 3 . Rising q, while keeping d\ and di con- 
stant, allows decreasing N s while maintaining the same gap 
properties. We illustrate this result for a single satellite with a 
series of three calculations, referred to as Runs a, b, and c. In 
all three models, d\ = 8 and di — 100. The main parameters 
of these runs are summarized in Table Q] The disk aspect ra- 
tio increases from h = 1% (Run a) to h = 5% (Run c), which 
comes to increase the satellite to primary mass ratio q and the 
kinematic viscosity v by about two orders of magnitude. All 
runs have 7V r = 448. The location of the grid's inner and outer 
edges, r m j n and r max , is adjusted such that Ar = r max — r m ; n = 
20/?h (5ina X /^H = 4.5%). The number of cells along the az- 
imuthal direction is taken as N s = 47ia/38 max , it is thus ss 5 
times smaller in Run c than in Run a. In all these runs, the 
satellite is held on a fixed circular orbit in order to focus on 
the gap properties at a given semi-major axis. The initial disk 
density at the satellite location is 10~ 4 in code units. Since 
the gas self-gravity is neglected, and the satellite is on a fixed 
orbit, the gas density value can be chosen arbitrarily. 

Azimuthally-averaged density profiles at 3000 orbits are 
plotted against A/Rn in the left panel of Figure [T] At a given 
radius, azimuthal averages exclude the disk parts located less 
than a few Hill radii away from the satellite. In scaled units 
A/Rn, the width and depth of the gaps are in good agreement. 
Although not shown here, we also checked that the density 
distribution inside the satellite's Hill radius is almost indistin- 
guishable from one calculation to another. This is an impor- 
tant point to keep in mind since, as will be shown below, the 
evolution of binary satellites is primarily driven by the density 
distribution in the satellite's Hill radius. 

We display in the right panel of FigureQ]the time evolution 
of the specific torque (torque per unit satellite mass) exerted 
by the disk on the satellite, for our series of three runs. The 
specific torques differ by a factor of order unity. The slight 
increase of the torque amplitude from Run a to Run c stems 
from the fact (i) the surface density at the outer horseshoe 
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Figure 1. Results of the calculations of §|5] whose main parameters are summed up in Table[T] Left: azimuthally-averaged density profiles at 3000 orbits. The 
azimuthal averaging excludes the disk parts located less than a few Hill radii away from the satellite. The x-axis displays the disk-satellite radial distance A, in 
units of the satellite's Hill radius. The initial surface density at the satellite's location equals 10~ 4 in code units. Right: time evolution of the specific torque 
exerted by the disk on the satellite. 



separatrix, namely at A ss 2Rn, takes very similar values in all 
runs, whereas (ii) the density of gas left in the gap, which 
notably contributes to the positiv e horseshoe drag, slightly 
decreases from Run a to Run c (Crida & Morbidellil I2007L 
their equation 15). We comment that the timescale to reach 
a steady-state varies between the simulations, which results 
from different gap's depletion timescales. A gap progres- 
sively deepens as fluid elements leave the satellite's horseshoe 
region. The timescale to build up a stationary gap profile is 
therefore proportional to the libration period of the horseshoe 
fluid elements, 



where x s is the half-width of the horseshoe region. For large 
values of the dimensiq nless parameter q/h?, x s ps 2Rh °< a 1 / 3 
(see e.g.,|Masset et al. 2006). The timescale to get a station- 
ary gap profile thus decreases with increasing q, and it should 
scale with h^ 1 . This is in good agreement with our results. 
The torques differ indeed by 5% from their final value at about 
2000 orbits for Run a, 1000 orbits for Run b, and 600 orbits 
for Run a. This comparison shows that, not only our rescaling 
method allows reducing the simulations computational cost 4 , 
it also helps to reach a faster steady-state, while getting very 
similar results in terms of the density distribution and specific 
torque. 

6. RESULTS OF CALCULATIONS OF THE RESCALED 
PROBLEM 

We have shown in §[5]that the original problem described in 
§|2] can be transformed into an equivalent, less computation- 
ally demanding rescaled problem, by using the two dimen- 
sionless parameters in the gap-opening criterion. By equiv- 
alent, we mean that the properties of the gap opened by the 
satellite (width, depth), and thus the value of the specific 
torque on the satellite are fairly preserved with the rescaling 
method. This result was shown for a single satellite. We will 
show in § l6.1l that the gap properties of single and binary satel- 
lites are in good agreement, allowing us to apply the rescaling 

4 We have checked that the computational time scales proportional to N s , 
as expected. 



method to binary satellites. All the calculation results pre- 
sented in this section therefore use the rescaling method. We 
come back to our original problem in §|7] 

Before presenting our results, we briefly enumerate the pa- 
rameters of the rescaled problem that we tackle in this sec- 
tion. The satellite to primary mass ratio is q = 10~ 3 , and 
the disk aspect ratio at the satellite location is h(a) = 5%. 
The disk's kinematic viscosity is constant, its value can be 
related to an equivalent alpha parameter at the satellite loca- 
tion a(a) = 4 x 10~ 3 . As shown in §|5] this set of parameters 
gives similar results (gap properties, specific torque on a sin- 
gle satellite) to that of our original, Galactic-center motivated 
problemin §0 q = 10~ 5 , h{a) < 1.1%, and a{a) = 10~ 3 . The 
initial surface density at the satellite location is 10~ 4 in code 
units. The initial Toomre parameter at the binary's loc ation 
~ 160, but (much) smaller values will be investigated in § 16.31 
The grid now includes N T = 800 rings and N s = 3000 sectors, 
and it extends from r m \ B = 0.55a to r max = 1.6a along the ra- 
dial direction. The satellite's Hill radius is resolved by about 
50 cells along the radial direction, and by about 33 cells along 
the azimuthal direction. Initially, the binary's semi-major axis 
is abin.o = 2 x 10~ 2 ao rs 0.3/?h- The softening length of the 
stars potential is e = abi n ,o/2, and, to avoid a violent initial 
evolution, the mass of the binary stars is slowly increased un- 
til it reaches its proper value at 10 orbital periods. 

Our calculations are performed in two steps. The satellite is 
first held on a fixed circular orbit for 500 orbits, during which 
the satellite opens a gap. During this stage, the force exerted 
by the disk on the satellite is ignored. Both the quantities 
a and at,i n therefore remain stationary during this first stage. 
Because initially flbin = 0.3/?h, the increase of the binary's 
eccentricity due to the central star remains small (on average, 
<?bin ~ 0.04). After 500 orbits, the force exerted by the disk on 
the satellite is included, inducing the satellite's orbital evolu- 
tion. For comparison, calculations were performed both with 
a binary and a single satellites of the same q. 

6.1. Binary satellite on a fixed circular orbit 

We describe in this paragraph the results obtained in the 
first 500 orbits the binary satellite is held on a fixed circular 
orbit. The disk surface density obtained at 500 orbits is dis- 
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Figure 2. Calculation result of the rescaled problem (q = 10 , h = 5%, a = 4 x 10 ). Top-left: surface density contour at 500 orbits of the gaseous disk 
surrounding the supermassive black hole, located at x = 0, y = 0. The binary's center of mass is located at x = 1, y = 0. Top-right: close-up of the density 
distribution around the binary (note that polar coordinates are used in this panel, the binary's center of mass is located at r = ft,, (p = <Pb). The dashed circle shows 
the binary's Hill radius. Part of the computational grid is overplotted in the bottom-left part of the panel. The bottom panels show a sequential time evolution of 
the disk density in the binary's Hill radius over almost half a revolution of the stars around their center of mass (which corresponds to about 0.05 orbital periods 
around the central black hole). The color scale and the axes are the same as in the top-right panel. 
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Figure 3. Left: surface density profile time-averaged over the 500 orbits the satellite is held on a fixed circular orbit. In contrast to the left panel of Figure [T] 
the azimuthal averaging takes into account the full azimuthal extent of the grid. The bottom x-axis displays the disk-satellite radial distance A, in units of the 
satellite's Hill radius. The top x-axis shows the separation from the central object. Right: running-time averaged specific torque exerted on each satellite (the two 
stars in the binary satellite case). Time is given in orbital periods of the satellite around the central black hole. 
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played in the top-left panel of Figure [2] The central object 
is located at x = 0, y = 0, the binary's center of mass is at 
x = 1, y = 0. A close-up on the density distribution around 
the binary is shown in the top-right panel of Figure [2] The 
stars are located at the center of the two densest points in this 
plot. The dashed circle depicts the satellite's Hill radius, the 
size of which is sw 1.4H(a). Part of the computational grid 
is displayed in the bottom-left part of this panel. The binary 
is prograde, stars thus rotate counterclockwise in this panel. 
We see that each star is lagged by a spiral wake (or tail). The 
impact of t hese wakes on the binary's orbital evolution is pre- 
sented in § 16.21 The bottom panels in Figure|2]display a time 
evolution sequence of the density inside the binary's Hill ra- 
dius during almost half a revolution of the stars around their 
center of mass (that is during m 0.05 orbital periods around 
the central black hole). The color scale and the axes are iden- 
tical as in the top-right panel of this figure. This sequence 
shows that the density in the trailing tails is strengthened as 
the stars cross the tidal wake induced by the binary's center of 
mass. This density enhancement arises from a larger velocity 
difference between the stars and the gas at the location of the 
tidal wake. 

Far enough from the binary's Hill radius, the disk density 
much resembles the typical disk density obtained with a sin- 
gle satellite, as expected (the quadrupole component of the 
binary's potential becomes negligibly small compared to the 
monopole outside of the binary's Hill radius). The left panel 
of Figure [3] compares the time-averaged surface density pro- 
files obtained with the single and binary satellites. Time- 
averaging is done over the 500 orbits of the simulations. In 
contrast to the left panel of Figure Q] the azimuthal averaging 
does not discard the material close to the satellite's location. 
This is meant to highlight that the disk mass in the vicinity 
of the single satellite is about two orders of magnitude larger 
than the one around the binary. In the binary case, the fast 
relative motion of the stars and of the background gas hin- 
ders gas accumulation inside the satellite's Hill radius. This 
is in contrast to the single satellite, where a mass ive circum- 
satell ite disk forms inside the Hill radius (see e.g., Cri da et all 
2009). Despite the large mass difference inside the Hill ra- 
dius, the width and depth of the gaps opened by the single 
and binary satellites are in good agreement. 

The right panel of Figure [3] displays the running-time av- 
eraged specific torque exerted on each star of the binary sys- 
tem (solid curves). Recall that time is in units of the satel- 
lite's orbital period around the black hole (as in all other fig- 
ures). Both torques are indistinguishable. For comparison, the 
running-time averaged specific torque on the single satellite is 
also shown (dashed curve). At 500 orbits, the torque obtained 
in the binary case is about 10% more positive, which results 
from a slightly less depleted gap, as illustrated in the left panel 
of Figure [3] In addition, we notice from the time evolution of 
the torques that the gap's depletion timescale takes somewhat 
longer with a binary satellite. In particular, the running-time 
averaged torque on the binary ha s no t reached a steady-state 
at 500 orbits. We will show in § 16.31 that it does not change 
our results for the binary's orbital evolution. 

6.2. Binary's orbital evolution 
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The simulation presented in § I6.1l was restarted at 500 or- 
bits, allowing the binary stars to feel the force exerted by the 
disk. As previously, we compare calculation results of single 
and binary satellites. The semi-major axis of the binary's cen- 
ter of mass (solid curve), and of the single satellite (dashed 




0.980 



0.020 



0.018 



o 
re" 
£ 

E 

</) 

in 

■>. 
i_ 
ra 
c 

m 



0.016 



0.014 



0.012 



0.010 



500 



520 540 
Time (orbits) 



560 




500 



520 540 
Time (orbits) 



560 




500 



520 540 
Time (orbits) 



560 



Figure 4. Orbital evolution of the binary star in its parent gaseous disk. The 
top panel shows the semi-major axis a of the binary's center of mass (solid 
curve). The calculation result with a single satellite is overplotted for com- 
parison (dashed curve). The binary's semi-major axis is displayed in the 
middle panel. The calculation was halted when the separation between the 
binary stars became approximately smaller than the stars softening length, 
depicted as an horizontal dash-dotted line in this panel. Recall that both a 
and Obi,, are expressed in units of ciq, the value of a before restart. The bi- 
nary's eccentricity e^ m is shown in the bottom panel. 
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Figure 5. Impact of the disk's density distribution inside the binary's Hill 
radius on its ha rden ing. We show the result of four runs restarted from the 
simulation of § 16. II In these four runs, the calculation of the force (or the 
torque) exerted by the disk on the binary stars excludes a fraction / of the 
binary's Hill radius Ru (centered on the binary's center of mass). The fraction 
/ equals 0, 1/4, 1/2 and 1 from bottom to top. The particular case with no 
material exclusion (/ = 0) corresponds to the calculation results shown in 
Figure [4] When the torque exerted by all the fluid elements in the binary's 
Hill radius is discarded (f = 1), the binary does not harden with time. The 
semi-major axis of the binary's center of mass is displayed in the inset plot 
for all four calculations. 

curve) are displayed in the top panel of Figure |4] Both the 
single and binary satellites migrate inward. The binary's mi- 
gration rate is slightly smaller, however, and it tends to de- 
crease with time. This behavior is consistent with the slow 
decrease of the running-time averaged specific torque shown 
in the right panel of Figure [3] and is due to the minor differ- 
ences in the gas density profiles depicted in the left panel of 
Figure|3] Extrapolating the results in the top panel of Figure|4j 
we find that the migration timescale of the binary's center of 
mass is about 5000 orbits, or ~ 8 x 10 6 yrs at 0.1 pc. This 
is about a factor of 2 shorter than the prediction of Eq. ©, 
which however estimates the minimum timescale expected in 
the type II migration regime. We have checked that this dis- 
crepancy arises from the torque exerted by the circum-satellite 
disk inside of the satellite's Hill radius. Another restart run 
with the single satellite, wherein the calculation of the force 
exerted on the satellite excludes the disk material inside of the 
satellite's Hill radius, gives a migration rate that is about half 
that of the simulation shown in the top panel of Figure H] and 
which is thus consistent with the pace expected in the disk- 
dominated type II migration regime. This behavior is further 
commented and illustrated for the binary satellite in the inset 
panel of Figure [5] 

As shown in the middle panel of Figure [4] the binary's 
semi-major axis abin decreases with time. We observe that the 
amplitude of the hardening rate, namely the quantity |flbin|. 
slowly increases with time. The calculation was stopped when 
the separation between both stars went slightly below the soft- 
ening length £ = abin,o/2, depicted as a dash-dotted line in 
this panel. Below this separation, the binary satellite mostly 
interacts with the surrounding disk like a single-star satellite. 
In particular, we have checked that a massive, circular cir- 
cumbinary disk starts to form around both stars, like with a 
single-star satellite. Our calculation result indicates that the 
binary's semi-major axis is reduced by a factor of 2 in m 70 
orbits around the central black hole, or in about 1000 orbits 
of the stars around their center of mass. This corresponds to 



~ 10 5 yrs at 0.1 pc in the Galactic center case. In this partic- 
ular simulation, the hardening timescale of the binary is ~ 40 
times shorter than its migration timescale. 

The bottom panel of Figure |4] shows the time evolution of 
the binary's eccentricity ebin- The interaction with the gas 
slowly decreases both the averaged value, and the amplitude 
of the eccentricity's oscillations. Although not shown here, 
we have checked that the eccentricity of the binary's center of 
mass remains negligibly small with time. 

The binary's hardening results from the formation of a spi- 
ral tail at the trailing side of each star (lEscala et al.l 120041 
iKim et al.ll2008b IStahler!l2010f) . iKim et all d2008l) showed that 
each perturber is dragged backward by its own wake, and 
pulled forward by the wake of its companion. Differently 
said, each companion of a binary system loses angular mo- 
mentum due to its own wake, and gains angular momentum 
from its companion's wake. For equal-mass perturbers, the 
ratio of these positive and negative torques is always smaller 
than unity, and it varies with t he perturbers' v elocity with re- 
spect to the background gas jKim et al.l l2008). The net effect 
of the wakes is therefore to extract angular momentum from 
each star, which is ultimately transmitted to the disk though 
viscous shear or shocks. This is analogous to the net effect of 
the inner and outer wakes generated by the binary's center of 
mass in the disk around the central black hole. 

To further illustrate the connection between the hardening 
of the binary and the flow inside its Hill radius, we performed 
additional simulations wherein the calculation of the force ex- 
erted by the disk on the stars excludes a fraction / of the bi- 
nary's Hill radius. Differently said, the binary stars feel the 
force exerted by all the fluid elements, except those located 
inside a circle centered on the binary's center of mass and of 
radius fRa- We took three values for /: 1 /4, 1 /2 and 1. The 
results of these simulations are displayed in Figure [5] The 
result with / = 0, which corresponds to the simulation pre- 
sented in Figure@] is overplotted for comparison. The time in 
x— axis is that after the restart (namely, after the first 500 fixed 
circular orbits during which the satellite opens a gap). The bi- 
nary's hardening rate decreases with increasing /. When the 
torque evaluation discards all of the Hill radius, the binary's 
hardening rate is vanishingly small. This result clearly shows 
that the hardening of the binary is caused by the gas inside its 
Hill radius. In addition, the very similar results obtained with 
/ = and f =1/4 show that the hardening does not stem from 
the gas directly surrounding the stars (the two densest points 
in the top-right panel of Figure |2), which is poorly resolved. 
The progressive decrease of the hardening rate from / = 1/4 
to / = 1 confirms that the hardening does arise from the tails 
lagging the stars. 

The inset plot in Figure [5] depicts the time variation of the 
semi-major axis of the binary's center of mass obtained with 
the previous simulations. The migration rates obtained with 
/ = and / = 1/2 are in very close agreement, whereas the 
binary's hardening rates differ by a factor of ~ 2. It suggests 
that the hardening of the binary has a very limited impact on 
its migration rate. For / = 1 (no hardening), the migration rate 
takes a smaller valu e, in agreement with the findings with a 
single-star satellite (Crida et al. 2009). Our results thus show 
that the hardening rate of the binary is primarily triggered by 
the disk inside the binary's Hill radius, whereas its migration 
rate is controlled by the disk outside of its Hill radius. Our cal- 
culations do not suggest there is a significant feedback from 
one rate on the other. This result does not mean that mod- 
elling the global disk structure around the central black hole 
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is unnecessary. Akin to the single satellite case, the flow of 
gas entering and leaving the Hill radius depends indeed on 
the interaction of the satellite with the whole disk (e.g., by 
the capture of fluid elements inside the satellite's horseshoe 
region). We finally comment that as the binary hardens, more 
gas can flow inside the binary's Hill radius, which can explain 
the slow increase of the binary's hardening rates seen in the 
middle panel of Figure |4] 

6.3. Parameter-space study 

We have shown in § 16.21 that the interaction between a bi- 
nary star and the gaseous disk it is embedded in leads to the 
hardening of the binary as it migrates inward. The harden- 
ing timescale is found to be much shorter than the migration 
timescale. This section is aimed at investigating the depen- 
dency of the hardening timescale upon the disk and satellite 
properties, and upon numerics. A detailed exploration of the 
parameters space is outside the scope of this paper. Disk and 
satellite parameters are allowed to vary within a narrow range 
of values such that the main assumption of our model (the 
binary opens a gap) is fulfilled. 

We first consider the impact of the unperturbed surface den- 
sity at the binary's location. Two additional simulations were 
performed with £o( fl o) = 5 x 10~ 4 and Eo(ao) = 2 x 10~ 3 . 
The initial Toomre parameter go at the binary equals 32 and 
8, respecti vely. The result of these simulations, and that of 
the run of § 16.21 are displayed in the top-left panel of Figure [6] 
(solid curves). Increasing the unperturbed density results in 
a faster hardening of the binary. For (2o( fl o) = 8, which can 
be seen as a representative value for our original model in 
§ |2] the binary's semi-major axis is reduced by a factor of 2 
in only ~ 5 orbits. In that case, the mass enclosed in the bi- 
nary's Hill radius when it starts to harden is about 50 times 
smaller than the binary's mass itself. The dependence of the 
hardening rate with varying Eo(ao) can be interpreted as fol- 
lows. As the binary is held on a fixed orbit, and because 
self-gravity is neglected, the fractional change (E — Eo)/Eo 
of the gas surface density is independent of Eq. This is true 
in particular inside the binary's Hill radius, which suggests 
that the net torque exerted by the trailing tails on each star, 
and therefore the binary's hardening rate |dbin|i should scale 
with Eo(ao). The upper dashed curve displays a linear fit to 
«bin( f ) f°i" £o( fl o) = 10~ 4 . The two bottom dashed curves de- 
pict previous fit with a slope increased proportionally to the 
unperturbed density. The good agreement with the calcula- 
tion results confirms that |dbin| is approximately proportional 
toEo(ao)- 

We have also checked the dependence of the hardening rate 
o n the restart time of our simulations. Restarting the run of 
§ 16.11 at 280 orbits gives the result overplotted in the top- 
left panel of Figure [6] as a dotted curve . It is in very close 
agreement with the calculation result of § l6.2l where the restart 
time is 500 orbits, despite the fact the running-time averaged 
torque has not yet reached a steady-state in either case (see 
right panel of Figure|3}. This agreement is consistent with the 
fact that the binary's hardening is driven by the density distri- 
bution within its Hill radius, which attains a faster steady-state 
than the density distribution in the gap. We have obtained the 
same good agreement with Eo(flo) = 2 x 10~ 3 . 

We now assess the impact of the parameters altering the gap 
properties, taking Eo(flo) = 2x 10~ 3 for convenience. The re- 
sult of calculations with varying the disk's kinematic viscos- 
ity is displayed in the top-right panel of Figure|6](the run with 



V = 10~ 5 , which corresponds to a = 4 x 10~ 3 at the planet's 
location, is the same as in the left panel of this figure). In- 
creasing the viscosity results in a faster hardening. Again, 
this result can be explained by the fact a larger viscosity in- 
duces a larger mass flow in the satellite's Hill radius, thereby 
increasing the amplitude of the drag force on each star. As in 
the top-left panel of Figure|6] the dashed curves highlight that 
the binary's hardening rate is proportional to the disk viscos- 
ity. We have checked that this result is quantitatively consis- 
tent with the density in the binary's Hill radius scaling propor- 
tional to the viscosity. We have also checked that, as expected, 
the migration rates obtained in these two runs scale with the 
disk viscosity. 

The bottom-left panel of Figure [6] compares the results ob- 
tained with a disk aspect ratio /; = 5% and h = 4%. Our 
simulations show that the hardening is slightly faster with a 
smaller disk aspect ratio, the hardening timescales differing 
by about 20%. In the same vein, the bottom-right panel of 
Figure [6] shows results with different (albeit close) values of 
the satellite to primary mass ratio, ranging from ^ = 6x 10~ 4 
to q = 1.5 x 10~ 3 . These simulations have the same abi n ,o- 
Our calculation results show that, for the range of values that 
we considered, q also has little impact on the binary's hard- 
ening rate. We have checked that the density distributions in 
the binary's Hill radius show very little dependence with the 
aforementioned values of h and q, which is consistent with 
our findings. The gap's depth, however, largely differs be- 
tween these runs, and so do the migration rates. We point 
out that in the simulation with q = 6 x 10~ 4 and h = 5%, the 
bina ry experiences inward runaw ay migration, in agreement 
with lMasset & Papaloizoul (120031) . 

We now discuss the dependence of our findings o n the res- 
olution of our calculations. As already stressed in § 14.21 res- 
olution is controlled by the softening length £ of the stars po- 
tential, and by the grid's resolution. We study in this para- 
graph the impact of both parameters. Two additional calcula- 
tions were performed with Eo(flo) = 2x 10~ 3 : one with dou- 
bling the number of grid cells along each direction, and one 
at our standard resolution but with e = 0. 1H (that is, half our 
standard softening length). Results are depicted in Figure [7] 
For our fiducial value of e, increasing the grid's resolution is 
found to have no impact on the hardening rate. Decreasing 
the softening length at same resolution yields however a dra- 
matic raise of the hardening rate (dotted curve). While the 
density inside the binary's Hill radius, in particular at the lo- 
cation of the tails, takes similar values for both values of £, 
the gas density at the immediate vicinity of the stars is in- 
creased by about one order of magnitude for e = Q.IH. It 
clearly suggests that the large mass accumulation around the 
stars is severely under-resolved, in contrast to the results pre- 
sented with our fiducial softening length, and that the resulting 
increase of the hardening rate is a numerical artifact triggered 
by a lack of resolution. To illustrate this statement, we per- 
formed an additional simulation at e = 0.1H, wherein the bi- 
nary's time-evolution discards the torque exerted by the fluid 
elements located inside a circle of radius Rn/4 ce ntere d on 
the binary's center of mass. As already stated in § 16.21 such 
radius is meant to discard the torque due to the mass accumu- 
lation around the stars, with very small impact on the torque 
exerted by the tails (see Figure |5j. The result of this simula- 
tion is overplotted in the bottom-right panel of Figure [7] and 
shows a good agreement with the result at larger softening. 
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Figure 6. Influence of the disk and satellite parameters on the time evolution of the binary's semi-major axis. In all panels, the horizontal dash-dotted line shows 
the value of the stars softening length. Top-left: results of simulations with increasing the unperturbed surface density £o(<?o) at the binary's initial location (solid 
curves). The black solid curve shows the result of the simulation in 5 I6.2I The black dotted curve shows the result of the same simulation restarted at 280 orbits 
(instead of 500 orbits). Dashed lines show linear fits to the numerical results, assuming the hardening rate |dbinl at the restart time scales with Eot^o) (see text). 
Top-right: results of an additional simulation with a smaller viscosity. Again, dashed lines show linear fits to the numerical results, assuming |dbinl scales with v. 
Bottom-left: result of a simulation with a slightly smaller disk aspect ratio. Bottom-right: calculation results with different satellite to primary mass ratios. 
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Figure 7. Dependence of the binary's hardening rate with numerics. We 
show the impact of the grid's resolution at a fixed softening length (top solid 
curves), and of the softening length at a fixed grid's resolution. In the latter 
case, results are shown with discarding the torque of the fluid elements lo- 
cated inside 1/4 of the binary's Hill radius (bottom solid curve), and without 
such exclusion (dotted curve). The corresponding values of e are depicted as 
horizontal dash-dotted lines. 



6.4. Case of a retrograde binary 

The calculation results presented in previous sections as- 
sume a prograde binary, with the angular momentum vectors 
of the disk and of the binary being aligned. It is possible that 
this relative inclination takes very different values when the 
binary forms in the disk, or after it formed, through dynami- 
cal interactions with other stars or more massive objects. An- 
other possibility is that the binary did not form in the disk, 
and was captured. One expects that the gaseous disk damp- 
ens the relative inclination on a rather short timescale. Since 
our calculation results suggest that the binary's hardening may 
also occur on short timescales, it is relevant to investigate the 
general case of a random disk-binary relative inclination. We 
limit ourselves to the case of a retrograde binary. The binary's 
center of mass still has a prograde motion around the black 
hole. 

We performed one simulation with a retrograde binary, us- 
ing the procedure and parameters described at the beginning 
°f §0 Akin to the prograde case shown in Figure|2j the disk's 
surface density obtained at 500 orbits with the retrograde bi- 
nary is displayed in the left panel of Figure [8] (note that the 
color scales differ in both figures). We observe that the gap 
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Figure 8. Left: surface density contour at 500 orbits of the disk perturbed by a retrograde binary held on a fixed circular orbit. Right: close-up of the density 
distribution around the retrograde binary. As in Figure|2] the dashed circle shows the binary's Hill radius, and part of the computational grid is overplotted in the 
bottom-left part of the panel. 
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Figure 9. Orbital evolution of retrograde and prograde binaries. Hardening rates and migration rates are compared in the left and right panels, respectively. The 
horizontal dash-dotted line in the left panel shows the value of the stars softening length. 



is less deep in the retrograde case. We have measured that 
the azimuthally-averaged density in the gap is a factor of 2 to 
3 larger compared to the prograde case. The right panel of 
Figure [8] shows a close-up of the disk density around the ret- 
rograde binary, whose Hill radius is drawn as a dashed circle. 
In this panel, the stars rotate clockwise. A comparison with 
the top-right panel of Figure |2]reveals that the mass enclosed 
in the Hill radius is about one order of magnitude larger in 
the retrograde case. We attribute these differences to a larger 
velocity difference in the retrograde case between the stars 
and the fluid elements approaching the stars on horseshoe tra- 
jectories. Those fluid elements are therefore less deflected as 
they get closer to the stars, and they can enter the binary's Hill 
radius more easily. 

Similarly as in the prograde case, overdense spiral tails lag 



the stars. The tails should therefore extract angular momen- 
tum from the retrograde binary, which we expect to harden 
with time, just like the prograde binary. This is confirmed in 
the left panel of Figure [9] where the time evolution of the bi- 
nary's semi-major axis is displayed for both the prograde and 
the retrograde binaries. The horizontal dash-dotted curve de- 
picts the stars softening length. The hardening timescale is 
a factor of 5 to 10 shorter than in the prograde case, which 
can be attributed to the relative difference of the density in the 
Hill radius between the retrograde and prograde binaries. The 
right panel of Figure [9] compares the semi-major axis of the 
binary's center of mass in both cases. The inward migration 
of the retrograde binary is much slower on average, which is 
most presumably due to a shallower gap profile. 
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Figure 10. Surface density contour at 300 orbits obtained with the parame- 
ters of the original problem (q = 10 -5 , h=l%, and a = 10~ 3 ). The binary's 
Hill radius is shown as a dashed circle, and part of the computational grid is 
overplotted in the bottom-left part of the panel. 
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Figure 11. Time evolution of the binary's semi-major axis obtained with the 
disk and satellite parameters of the original problem (solid curve), and of the 
rescaled problem with same initial density (dashed curve). For comparison 
purposes, the ratio t/bin /«bin,o is displayed in y-axis. The time in x-axis is in 
units of the orbital period of the binary's center of mass around the central 
black hole. The horizontal dash-dotted line depicts the ratio e / c bin ,o> which 
is identical in both simulations. 



7. BACK TO THE ORIGINAL PROBLEM 

Our original problem is to investigate the tidal interaction of 
a massive (q = 10~ 5 ) binary star and its natal thin (h ~ 1%) 
gaseous disk. We have justified in § [5] that, for a single-star 
satellite, one can instead tackle a less computationally de- 
manding rescaled problem, where the properties of the disk- 
satellite interaction (gap properties, amount of gas in the satel- 
lite's Hill radius, value of the specific torque on the satellite) 
remain essentially unchanged. The calculation results in § [6] 
have been obtained with the parameters of the rescaled prob- 
lem: q — 10~ 3 , h — 5% and a = 4x 10~ 3 . We have shown 



that, in addition to migrating inward, the binary system hard- 
ens, due to the net drag force exerted by each companion's 
wake inside the binary's Hill radius. Our results indicate that 
the hardening timescale is inversely proportional to the unper- 
turbed density in the binary's Hill radius, and that it is shorter 
than the migration timescale by typically one to two orders of 
magnitude. 

We have shown in § 16. 1 I that the gap properties and the spe- 
cific torque obtained with single and binary satellites are in 
good agreement. This agreement implies that the steady-state 
migration rate of a binary satellite should be similar in the 
original and rescaled problems, as in the case of a single satel- 
lite detailed in § [5] However, it is not clear whether the hard- 
ening rate remains unchanged when applying the rescaling 
method. To check this, we performed a simulation with the 
parameters of the original problem. The unperturbed surface 
den sity at the binary's location is 10~ 4 , as in the simulation 
of § 16.21 This sets the initial Toomre parameter at the binary's 
location to 32. For comparison purposes, we kept the same 
ratios flbin/^H an d fl bin/£> which leads to a^ n « 4 x 10~ 3 a and 
£«2x 10~ 3 a. The satellite's Hill radius is resolved by about 
50 cells along the radial direction, and 24 along the azimuthal 
direction (to be compared with 50 and 33, respectively, in the 
simulations of § |6). A close-up of the gas surface density 
around the binary star is displayed at 300 orbits in Figure [lOl 
It is very similar to the density obtained at 500 orbits in the 
rescaled problem (top-right panel of Figure [2]). This simili- 
tude is due to (i) the use of (almost) the same dimensionless 
parameters in the gap-opening criterion, and to (ii) the fact the 
density in the satellite's Hill radius reaches a faster steady- 
state than the density in the gap. We notice that the density at 
each perturber's wake is a factor of w 2 larger in the original 
problem. 

The time evolution of the binary's semi-major axis obtained 
with previous simulation restarted at 300 orbits is depicted in 
Figure [TT] (solid curve). To compare with the result of the 
rescaled problem (dashed curve), the ratio abin/«bin,o is dis- 
played in y-axis. The binary's semi-major axis in the orig- 
inal problem decreases by a factor of two in about 8 to 10 
orbits, or in about 130 orbits of the stars around their center 
of mass. This corresponds to ~ 1.5 x 10 4 yrs at 0.1 pc in 
the Galactic center case. The slowing down of the binary's 
hardening rate after « 8 orbits is most presumably due to the 
fact that «bin becomes comparable to the softening length of 
the stars potential. The hardening timescale in the original 
problem is shorter than that of the equivalent rescaled prob- 
lem by a factor of 8 to 10. This timescale difference can be 
qualitatively interpreted as due to different mass ratios of the 
binary and of the gas in the binary's Hill radius. Between 
the original and rescaled problems, the binary's mass is in- 
creased by a factor of 100, whereas the mass inside the bi- 
nary's Hill radius in increased by a factor of ~ 20 (ratio of the 
Hill radii squared). This simple argument suggests a harden- 
ing timescale ~ 5 times larger in the rescaled problem, which 
is in correct agreement with our findings. 

Some furt her insight can be gained by using the analytic 
prediction of Stahler (2010). Through a linear analysis of the 
Euler equations, he derived the timescale for the orbital de- 
cay of a low-mass binary embedded in a three-dimensional 
static medium of gas. Applying his equation 56-a to our disk 
model (writing the sound speed c s = HQ., the 3D unperturbed 
gas density po ~ Ho/ 2H, noticing from our simulations that 
a typical density in the binary's Hill radius is the unperturbed 
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density at the binary's initial location), the extrapolated hard- 
ening timescale Th ar d may be estimated as 
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where T or \, is the orbital period of the binary's center of mass 
around the central black hole. Since our original and rescaled 
simulations have the same Q/h (that is, the same unperturbed 
density at the binary's location), Eq. dTTb predicts that the 
hardening timescale in the original problem should be a fac- 
tor of 8 snorter than in the rescaled problem, which is in good 
agreement with our results of simulations. For convenience, 
Eq. ( fTTT i can be recast as 
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where Mbi n denotes the binary's mass, and where T and £ are 
the unperturbed gas temperature and surface density at the bi- 
nary's location. In Eq. (fT2l . we have assumed a mean molec- 
ular weight n = 2.4, and an adiabatic index 7=5/3. 

With the parameters of the original problem, Eq. ( fTTT > pre- 
dicts a hardening timescale < 0.1 orbit, that is about two or- 
ders of magnitude shorter than in our simulation. We com- 
ment that, while Eq. (fTTT i may capture the right functional 
dependence of the hardening timescale with the disk and 
satellite parameters in our problem, a quantitative compari- 
son with our results of si mulations seem s inappropriate, since 
the key assumption in Stahler] d2010l) 's analysis (each per- 
turber's wake induces a small perturbation of the background 
gas density, which can be described with a linear analysis) is 
clearly not satisfied in our simulations (see e.g., Figure ITOb. 
Another source of discrepancy may also be triggered by the 
time-varying velocity difference in our simulations between 
the binary and the gas inside the Hill radius. 

8. DISCUSSION 
8.1. Motivation 

The primary astronomical motivation of the investigation 
presented here is to explore the origin of the "S-stars", a pop- 
ulation of relatively massive stars in the inner 0.05 pc of the 
Galactic center. The S-stars could be related to those in copla- 
nar orbits just outside this inner region. It is natural to assume 
that these d isk-stars may have eith er coherently formed in 
(lLevinl2007l) or rapidly accreted gas (lArtymowicz et al.l 19931) 
from a common gaseous disk. 

Most of the S-s tars have large eccentricities and random 
inclinations (e.g., iLu et al.l 120091) . In contrast, the nearby 
disk-stars have smaller eccentricities and their velocity dis- 
persion is an order of magnitude smaller than the local Kep- 
lerian speed. One potential explanation for their large eccen- 
tricities and random inclinations is that the S-stars originated 
from the tidal disruption of binary systems that were scat- 
tered to the proximity of t he central supermassive black hole 
on nearly parabolic orbits (Gould & Ouillen 2003). This sce- 
nario would also account for the observed hypervelocity stars 
in the Galactic halo, provided some of the scattered binaries 



have very close separat i on dHi lls 1988; Yu & Tremainell2003t 
iGualandris etail 12001 ILu et all 120 lOt iZhang et alJ 1201 Oh . 
typically within a few solar radii. 

In this paper, we consider a unified scenario for the young, 
massive stars near the Galactic center (both disk- and S-stars) 
and high-velocity stars in the Galactic halo. We assume their 
progenitors were binary stars that formed at large distance (a 
fraction of a pc) from the SMBH, in a gaseous disk that we 
expect existed in the Galactic center in the past. In princi- 
ple, similar to the solar neighborhood, binary systems could 
be as comm on as single stars and have a logarithmic period 
distribution (Duquennoy & Mayor 1991). Around the Galac- 
tic center, however, tidal interaction between binary systems 
and their natal disk induces angular-momentum exchange be- 
tween them. The main outstanding issue is whether these bi- 
naries' center of mass would migrate (presumably inward) ex- 
tensively before they harden substantially. 

8.2. Model 

We have investigated the tidal interaction of a binary star 
and its natal gaseous disk by means of two-dimensional hy- 
drodynamical simulations. For illustration purpose, we have 
considered a fiducial binary system, comprising two equal- 
mass (15 Mq) stars embedded in a thin (h = 1%), moder- 
ately turbulent (a = 10~ 3 ) gaseous disk around a 3 x 10 6 Mq 
SMBH. These initial and boundary conditions are such that 
the binary opens a deep gap around its orbit, and is thus sub- 
ject to type II migration. We have shown in § [5] that this gap- 
opening process has helped us transform our original problem 
into an equivalent, less computationally demanding rescaled 
problem. We have presented results of calculations of the 
rescaled problem in § [6] and of the original problem in § [7] 

We have made a number of simplifying assumptions in the 
models presented here (see §[2j. We have neglected the pro- 
cess of gas accretion onto the stars, assuming the latter al- 
ready reached their final mass at the beginning of the simu- 
lations. We have also assumed throughout this paper that the 
binary stars formed in a circular orbit, while the fragmenta- 
tion of an eccentric g aseous disk would put them in initially 
eccentr ic orbits (e.g., Nayakshin et al. 2007: [Alexander et alJ 
2008b). While we expect many stars to form at the same time 
in a fragmenting gas disk, we have considered the evolution 
of a single isolated satellite (one binary star). The presence of 
many embedded stars, each propagating density waves, and 
opening a gap for the most massive ones, would dramatically 
impact the gas disk structure. The properties of disk-satellite 
interactions, and the resulting migration timescale, could be 
very different from those of a single isolated satellite. The 
interaction between different satellites could excite their or- 
bital eccentricity and alter their migration rate, leading for ex- 
ample to ca ptures into mean-motion resonance or scatterin g 
events (e.g.. ICresswell & Nelson! 12008: Mar zari et al.ll2010h . 
Also, encounters with other single stars may help to further 
harden binary stars, especially at short sepa rations from the 
supermassive black hole (e.g., Perets 2009). We defer this 
interesting regime for future study. 

The properties of the gaseous disk after star formation are 
largely uncertain and poorly constrained. We have grossly 
simplified the evolution of the underlying gaseous disk, in par- 
ticular by neglecting the depletion of the gaseous disk, and the 
increase of its temperature under star formation. Local heat- 
ing of the gas disk by the intrinsic luminosity of the embedded 
stars could modify the mass flow inside the binary's Hill ra- 
dius, and therefore the binary's hardening rate. The progres- 
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sive depletion of the gas disk also implies a finite timescale 
over which disk-satellite interactions are e ffective. In the 
contex t of Sgr A*, we mention the work of Navakshi n et alj 
(2007), who performed hydrodynamical simulations of star 
formation in a gaseous disk around the central SMBH, includ- 
ing a crude treatment of thermal accretion feedback (a fraction 
of the star accretion luminosity is converted into disk heat- 
ing). Their simulations indicate that the disk's mass may be 
reduced by a factor of two in at most 10 5 yrs (about 60 orbital 
periods at 0.1 pc). It is much shorter than the timescale for 
gas accretion ont o Sgr A* through turbulent transport of an- 
gular momentum (Navakshin & Cuadra 2005). It is also much 
shorter than the expected timescales for type I and type II mi- 
gration, but not necessarily short compared to the timescale 
for type III migration (see § [3]). Still, this potentially short 
lifetime of the gas disk is a challenge for migration driven by 
disk torques. However, it may not be a stringent constraint 
on the hardening of massive binary stars driven by the gas 
disk, which we find to occur on a timescale comparable to, or 
s horte r than the aforementioned disk lifetime, as is recalled in 

§E3 

8.3. Summary of our results 

Our main result is that the binary gets harder as it migrates 
toward the central black hole. The hardening timescale is 
much shorter than the migration timescale, which is very sim- 
ilar to that of a single satellite of the same mass. The binary's 
hardening rate is essentially determined by the gas density in- 
side the binary's Hill radius. In that region, overdense spiral 
tails lag the stars, exerting a torque on them and hardening 
their orbit. The gas density at the tails location, and therefore 
the hardening rate, scale with the viscosity and with the un- 
perturbed surface density of gas at the binary's location (i.e. 
prior to gap formation). Assuming that at this location, the a 
viscous parameter is ~ 10~ 3 , and the initial Toomre parame- 
ter (before gap formation) is m 32, we find that the separation 
between the stars is reduced by a factor of 2 in about 10 orbits 
of the binary around the central black hole (see §|7]). Scaled 
to the Galactic center case, that is in only ~2x 10 4 yrs at 
0.1 pc, which is several orders of magnitude shorter than the 
expected migration timescale. A larger initial surface density 
would yield an even shorter hardening timescale. In the rest of 
this section, we discuss the evolution of low-mass binary stars 
(namely, what could have happened prior to our initial condi- 
tion), and the subsequent evolution of massive binary stars 
(that is, what could occur after the end of our simulations). 

8.4. Hardening of low-mass binary stars 

It is likely that disk-satellite interactions significantly im- 
pact the binary's semi-major axis from the early stages of 
the binary's formation, well before the binary has attained 
its asymptotic mass and potentially opened a gap around its 
orbit. In that case, the wake of each star would induce a 
low-amplitude density perturbation in the background gas sur- 
rounding the stars. Th is situation would resemble the one in- 
vestigated by Stahlei] (120101) . The hardening timescale ex- 
trapolated from his analytic study to our disk model is given 
at Eqs. ( fTTT i and ( fT2l . Assuming the unperturbed gas disk 
has Q ~ 32, h ~ 1% (that is, E - 6 g cm~ 2 and T ~ 2000 
K at the binary's location), and taking a low-mass binary 
(1 Mq stars, so that q/h 3 < 1 and the linearization assump- 
tion should be valid) with a^ ~ 0.37?h ( as m our simulations; 
here it corresponds to ctbin ~ 40 AU), Eqs. (fTTT i or (fT2l give 



Thard ~ 1 -2 x 10 4 yr at a = 0.1 pc. This timescale is similar 
to the one obtained in our 2D simulations, but note the strong 
dependence of Thard with both q and h. The dynamical evo- 
lution of a low-mass binary subject to hardening during its 
type I migration deserves follow-up investigation with self- 
consistent high-resolution 3D simulations. 

8.5. Subsequent binary evolution 

The results in §|6]and|2]provide the initial hardening rate for 
massive (gap-opening) binary stars with separation close to 
their Hill radius. We underscore that our simulations cannot, 
and do not address the final outcome of the binary's harden- 
ing, i.e. whether binary stars will ultimately coalesce, form a 
very hard binary, or become disrupted. We can only extrapo- 
late some possible implications of the binary's fast hardening 
in our model. 

8.5.1. Open ing of a cavity 

As a massive binary star shrinks, residual gas is preserved 
in a viscous disk between its orbit and its Hill radius. Tidal 
interaction between the binary and the disk around it leads to 
angular momentum transfer from the stars' orbits to the inner 
disk region. Similarly, angular momentum is removed from 
the outer region of the disk by the SMBH's tidal torque. The 
binary's hardening rate is determined by an equilibrium (out- 
ward) flux of angular momentum. Our results suggest this 
rate may be maintained at a fairly large value at least un- 
til the binary separation is less than about one sixth of its 
Hill radius. Our results also suggest that the disk interior to 
the binary stars' orbit should be cleared by their tidal torque 
such that this region is expected to evolve into a cavity. In 
that case, the inside of the binary's Hill radius would com- 
prise the binary star, the cavity where the binary now evolves 
in, and the surrounding gas disk referred to as the circumbi- 
nary disk (the cavity should not be confused with the gap that 
the satellite opens or has opened in the global disk structure 
around the SMBH). It remains however to be clarified how 
the timescale associated with gas removal from inside the bi- 
nary's orbit compares with the binary's shrinking timescale 
before a cavity is opened. In the next paper of this series, 
we will present high-resolution hydrodynamical simulations 
to make this comparison. If the binary manages to main- 
tain a cavity, the tidal interaction between the binary and 
the circumbinary disk would keep on hardening the binary, 
at slower rates than the initial pace we found here. Other- 
wise, the binary could keep on hardening until it coalesces. 
It is thus possible that some of the massive stars near the 
Galactic center result from coalescences of binary-star sys- 
tems, which might contribute to account for the top-heavy ini- 
tial mass function ob served in the two s tellar disk structures in 
the Galactic center (Bart ko et al.ll20101) . However, such stars 
would still be embedded in a gaseous disk around the SMBH, 
and t herefore in nearly-circula r orbits. Resonant relaxation 
(e.g., iRauch & Tremaindll996t iPerets et al J [20091) . the pres- 
ence of sho rt-lived massive stars or of intermediate mass black 
holes (e.g. . I Yu et al.l 120071). or an initially-eccentric gaseous 
disk (e.g., ICuadra et al.1 120081) could explain their observed 
current eccentricities. 

The hardening process with a cavity inside the binary's or- 
bit is very similar to type II migration: the tidal torque of the 
binary competes with the viscous torque from the circumbi- 
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nary disk (see e.g., IPringlel 1 199 It lArtvmowicz et al.1 1199 It 
ISver & Clarke! 1995tlIvanov et alj! 999). This problem shares 
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a number of analogies with the hardening of bi nary black- 
holes surrounded by a circumbinary disk (e.g., Cuad ra et al.l 
2009, where the transport of angular momentum is due to 
the disk self-gravity). We take the separation r a between the 
binary's center of mass, and the location in the circumbi- 
nary disk where most of the binary's angular momentum is 
deposited (that is, near the inner edge of the circumbinary 
disk), equal to abin- Further assuming that the gas density 
E at that location is still as large as ~ 10 2 g cm~ 2 , the ra- 
tio 47rEag in /M bin - 5 x 10~ 2 < 1 for a M bin = 30 M bi- 
nary with separation abin ^100 AU. The asymptotic hard- 
ening timescale ThaxAas of the binary thus corresponds t o th e 
timescale for satellite-dominated type II migration (see £j |3.21 i. 
Using Eq. (0, we have 



^hard,i 



67TVE' 



(13) 



where E and v denote the surface density and kinematic vis- 
cosity at r Q . Further writing v = ac 2 £2 _1 , with c s and Q. the 
sound speed and angular velocity of the circumbinary disk, 
we obtain 



Thard,as ~ 9 Myr X 



a 



10 



-3 



10 2 g cm 



-l 



( T \ 1 / M bin 

V2000K/ ^30M C 

Cbin \ ~ 3 / 2 

100 AU 



3/2 



(14) 



where we have assumed a mean molecular weight fi =2.4, 
and an adiabatic index 7=5/3. The unknown values of a, T 
and E near the inner edge of the circumbinary disk make the 
value of Thard as largely uncertain. In order to compare with 
the hardening timescale inferred from our simulations, we ar- 
bitrarily assume that a, E and T take similar values as in our 
physical model in §|2] which however describes the properties 
of the star-forming gas disk around the SMBH. This choice 
is also motivated by the simplifying assumption that the cir- 
cumbinary disk is in a quasi-steady state. The quantity vE in 
Eq. ( fT3] l then becomes independent of the separation from the 
binary's center of mass, and it may be evaluated near the outer 
edge of the circumbinary disk (where a, c s and Q. should not 
dramatically differ from their value in the star-forming disk 
about the central SMBH, at the binary's location). Thus tak- 



ing T ~ 2000 K, E - 10 g cm~ z , a ~ 10" J , M bin = 3OM , and 
flbin ^100 AU, we get Thatd as ~ 90 Myr. This timescale is, as 
expected, much longer than the hardening timescale with no 
cavity inside of the binary's orbit. Even with considering a 
more conventional larger a, smaller values of E due to the 
progressive disk depletion should maintain Thard as to a value 
that presumably exceeds a few Myr. This is comparable to, 
if not larger than the lifetime of the most massive stars near 
the Galactic center, indicating that stellar evolution should be 
taken into account (see § |8.5.2| >. 

We finally point out that there is no indication in Eq. ( fl4l > 
as to why the binary's shrinkage should stop. We mention 
the possibility that the gas in the circumbinary disk may 
be progressively stripped due to stellar winds, or to multi- 
ple close encounters w ith other stars, as was observed by 
iNavakshin et al.l (120071) . If the shrinkage is sufficiently effi- 
cient to decrease the binary's separation to a few stellar radii, 
mass transfer due to Roche lobe overflow will control the bi- 
nary's subsequent evolution. If the donor star expands faster 
than its Roche lobe, or shrinks less rapidly than its Roche lobe 
for a prolonged time, mass transfer will be unstable, which 



may lead to a common-envelope phase and, perhaps, a merger. 

8.5.2. Including stellar evolution: the explosive disruption 
scenario 

Including stellar evolution, and the likely possibility of 
unequal-mass binaries, we now consider the outcome of bina- 
ries that have reached small enough separations. If the more 
massive (primary) star evolves into a red giant with a radius 
comparable to the binary separation, both stars undergo com- 
mon envelope evolution. The engulfed binary hardens again 
due to the drag force exerted by the gas on each star. The en- 
ergy lost by the binary system leads to the expansion of the 
gas envelope. This is reminiscent of the opening of a cavity 
by a binary embedded in a gaseous disk, discussed in § 18.5.11 
The successful ejection of the common envelope, and there- 
fore the binary's survival, depend on the binary separation be- 
fore the common envelo pe stage, and on the prim ary's mass 
(see e.g., the review by Taam & Sandquist 2000). If the bi- 
nary survives the common envelope evolution (i.e. if it does 
not coalesce), it then comprises an Helium-rich star and an 
OB-star companion. 

Many stars near the Galactic center lose mass through Wolf- 
Rayet winds. With a substantial mass loss, the primary may 
evolve into Wolf-Rayet or blue supergiant stars. Mass loss 
typically results in increasing the binary separation, but in our 
scenario this trend is offset by the disk presence. In this case, 
the binary separation can be reduced to several solar radii, 
and may avoid another common envelope evolution if the pri- 
mary explodes as a supernova. The outcome of the primary's 
explosion depends on the fraction of the primary's mass that 
is ejected. If they explode, typical blue supergiants (with a 
mass ~ 10 — 20 Mq) may lose most of their residual envelope 
during their supernova explosion. If the mass lost to the su- 
pernova ejecta is more than half of the binary's mass before 
explo sion, the compact binary would become unbound dHillsl 
1983), and the secondary could attain a recoil velocity that is 
a significant fraction of the primary's surface escape velocity. 

What velocities can we expect in the supernova disruption 
scenario? The simplest approach is to look at the properties 
of the two-body orbit of the secondary about the central black 
hole. Before disruption, we assume that the binary's center 
of mass is on a circular orbit with semi-major axis a, and the 
binary's internal separation is «bin (we neglect the binary's 
internal eccentricity e b in f° r simplicity). We also assume that 
the primary's supernova has no impact on the secondary: the 
mass and velocity of the secondary before and after disruption 
are taken the same. We will come back to this assumption 
below. In the frame centered on the black hole, the velocity of 
the secondary is the sum of the velocity of the binary's center 
of mass about the black hole (v smb h), and of the velocity of 
the secondary about the binary's center of mass (Vbi n ). The 
angle 9 between both velocity vectors is a random variable 
distributed uniformly between and 2n, since in our model 
both velocities are coplanar. Their amplitude is given by 

/ Ttyf \ 1/2 / X -1/2 

WM -i I M smhh \ f a 
v S mbh ~ 390 km s x 



and 



Vbin ~ 133 km s x 



3.5 x 10 6 M, 



Mi \ ( M b 



0.1 pc 
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Figure 12. Outcome of the supernova disruption of a binary star orbiting a M sm t,h = 3 .5 x 10 6 M© black hole. Before disruption, the primary's mass is Mi = 20M ., , 
and the secondary's mass is either Mi = 15M;, (left panel), orM2 = 4.U . (middle and right panels). We assume that the primary loses enough mass so that the 
binary becomes disrupted. Before disruption, the internal and orbital angular momentum vectors of the binary are taken aligned. We assess in this figure the 
properties of the two-body orbit of the secondary with the black hole (see text for details). In all contour plots, the binary's internal separation at,in is depicted 
in x-axis, and the separation a between the black hole and the binary's center of mass is in y-axis. The left panel shows the mean eccentricity of the bound 
secondaries. The middle panel displays the fraction of escapers, namely the secondaries that become unbound to the black hole. The escapers' mean velocity at 
infinity is shown in the right panel. 



respectively. In Eq. ( fl6l l. M\ and Mi denote the mass of the 
primary and of the secondary, and M D i n —M\+ Mi is the bi- 
nary's mass. In the following, we take M sm b n = 3.5x 10 6 M Q , 
andMi=2OM . 

We first consider the case where Mi = 15 M©, a mass simi- 
lar to that of the S-stars. We vary a from 4 x 10~ 3 pc to 1 pc, 
and flbin from 0.1 AU to 4 AU (both with a logarithmic spac- 
ing). The minimum value of a D in is approximately the sum 
of the two stars radii. Its maximum value is set as about half 
the smallest Hill radius of our binaries sample. This choice 
ensures that, before the primary goes supernova, all binaries 
are bound with modest internal eccentricity. The left panel 
of Figure [12] shows a contour plot of the mean eccentricity of 
the secondaries that remain bound to the central object, as a 
function of abin (x-axis) and a (y-axis). The average is done 
over the different values of the angle 9. Large eccentricities, 
compatible with those of the S-stars, can be obtained with the 
supernova disruption scenario. We point out that, for suffi- 
ciently large a (at a given a D inX the mean eccentricity of the 
bound secondaries decreases, with a local minimum of ~ 0.6, 
and increases again. The presence of this local minimum is 
due to secondaries put on retrograde, close to circular orbits 
about the black hole. It is straightforward to show that the 
minimum's location satisfies Vb; n = 2v smD h, which can be re- 
cast as a ~ 0.6 pc x (a D in/0.1 AU) for our set of parameters. 

We next consider the case where Mi = 4 M Q , a mass sim- 
ilar to that of the hypervelocity stars, and we now vary «bin 
from 0.07 AU to 4 AU. The middle panel of Figure [121 shows 
the fraction of escapers, these secondaries that are energetic 
enough to become unbound to the black hole. The fraction of 
escapers increases with increasing a and decreasing abin, as 
expected. For maximally hard binaries, the fraction becomes 
unity for a > 0.4 pc. 

The escapers' mean velocity at infinity, measured in the 
frame centered on the SMBH, is shown in the right panel of 
Figure[T2] The large velocities that we obtain 5 suggest that the 
most energetic escapers could account for some of the Run- 

5 We have checked that the inclusion of the gravitational potential of the 
Galactic bulge, disk and halo (see e.g., Xue et al. 2008) has a negligible im- 
pact on estimating the escapers velocity at infinity. 



away Galactic OB-stars, as originally proposed by iBlaauwl 
(1961), or for some of the hypervelocity stars in the Galac- 
tic halo. The latter have typical Galactic rest-frame velocities 
between 400 and 700 km s _1 (e.g., iBrown et af] 120061) . Our 
results indicate that the ejection of secondaries that may still 
be going faster than ~ 400 km s~ 1 in the Galactic halo is pos- 
sible. It however requires binaries that are nearly maximally 
hard, and primaries that are sufficiently massive (the maxi- 
mum escapers velocity approximately scales proportional to 
\fM[; a primary's mass before explosion of ~ 10M© would 
yield a maximum escaper's velocity ~ 330 km s _1 ). 

The suggestion that the hypervelocity stars could orig- 
inate from the supernova disruption of binaries close to 
the SMBH is in appar ent contradiction with the results of 
iGualandris et al.1 (l2005h . who performed population synthe- 
sis studies and found that no binary can produce a hyperve- 
locity star by itself. However, in our scenario, binaries can 
be kept at very short internal separations due to the action of 
the gaseous disk. In contrast, standard population synthesis 
studies do not model the angular momentum loss of binaries 
in presence of a background disk, and thus produce few very 
hard binaries. We stress that our study mentions the possi- 
bility, but does not estimate the probability to get escapers 
with velocities comparable to those of the runaway and hy- 
pervelocity stars through supernova-disrupted binaries. We 
finally comment th at the distribution of the hypervelo city stars 
on two thin disks dLu et alj|20 10; Zhan g etaLll2010l) may be 
accounted for by both the tidal and the supernova disruption 
scenarios. In both cases, we need compact binaries originat- 
ing from two disks near the Galactic center. While the tidal 
disruption mechanism requires binaries on highly elliptical 
orbits, the supernova disruption mechanism requires binaries 
with internal and orbital angular momentum vectors that are 
nearly aligned. 

Only mean quantities are displayed in Figure [12] To give 
an idea of their standard deviations, we show in Figure Qj] 
the mean and standard deviation of (i) the eccentricity of 
the bound secondaries (black cross symbols, see left axis), 
and of (ii) the escapers velocity at infinity (red filled circles, 
right axis) resulting from nearly maximally hard supernova- 
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Figure 13. Eccentricity of the bound secondaries (black cross symbols, left 
axis), and escapers' velocity at infinity (red filled circles, right axis) result- 
ing from the supernova disruption of nearly maximally hard binaries with 
aligned internal and orbital angular momentum vectors. Error bars indicate 
the standard deviation of each quantity. 




0.01 0.10 
Semi-major axis at disruption [pc] 

Figure 14. Eccentricity (mean and standard deviation shown as error bars, 
left axis) and inclination (standard deviation, right axis) of the bound secon- 
daries resulting from the supernova disruption of nearly maximally hard bi- 
naries with misaligned internal and orbital angular momentum vectors. Two 
cases are contemplated: one for which the internal angular momentum vector 
of the binary before disruption points towards the black hole (black cross and 
square symbols), and one where it is parallel to v sm bh (blue filled circles and 
star symbols). 



disrupted binaries. 

The results shown in Figures [12] and Q~3] assume that the 
internal and orbital angular momentum vectors of the binary 
before disruption are aligned, resulting in a change of the sec- 
ondaries' eccentricity, but not of their inclination. Our model 
assumes that binaries are formed and stay in the same plane 
as that of the gas disk. While a more general scenario, in- 
cluding possible misalignments and their distribution, is out 
of this paper's scope, we briefly underscore the possibility 
to get secondaries with large inclinations through the super- 
nova disruption scenario. This statement is illustrated in Fig- 
ure[l4l where we display the eccentricity (mean, and standard 
deviation depicted as error bars), and standard deviation of 
the inclination of bound secondaries arising from the disrup- 
tion of nearly maximally hard binaries at different separations 
from the SMBH {M\ = 20M &1 M 2 = 15M and a bin = 0.1 
AU). Only two extreme cases are shown, one where the in- 



ternal angular momentum vector of the binary initially points 
towards the SMBH (black cross and square symbols), and one 
where it is parallel to v smbb (blue filled circles and star sym- 
bols). Our results show that it is possible to obtain secon- 
daries with large inclinations for a broad range of eccentric- 
ities (including large ones) through the supernova disruption 
of very hard, circular binaries. In contrast, the tidal disruption 
of such binaries can only lead to stars with small inclinations, 
since the ratio of the binary's orbital and internal velocities 

~ (Msmbh/Mrin) 1 '' 3 ^> 1- A statistical investigation aimed at 
comparing the outcome of the supernova disruption scenario 
with the orbital properties of the known S-stars will be pre- 
sented elsewhere. 

We stress that the above results assume a very simple dis- 
ruption model. In particular, we have discarded the impact 
velocity imparted by the ejected shell of the primary. Its ex- 
pression is given bv lTauris & fakeiisl (fl998h . Still assuming 
for simplicity that the secondary's mass is not altered by the 
supernova's blast wave, the impact velocity Vi mp can be esti- 
mated as 

^^G^) (17) 



where M sb and v s h denote the mass and ejection speed of 
the ej ected shell. Taking v s h ~ 10 4 km s -1 (iTauris & Takensl 
1998), Mi = 2OM , M 2 = 4M , and M sh - 12 M (such that 
the binary loses about half of its initial mass and becomes 
unbound), we get Vi mp ~ 1.5 kms -1 (ctbin/1 AU) -2 , that is 
Vimp ~ 0.7v bin (a bin /0.07 AU) -3 / 2 . We have checked that, for 
our range values of a and a b j n , the inclusion of the above im- 
pact velocity has a mild impact on our results (the maximum 
escapers' mean velocity is increased by 20%, to ~ 600 km 
s -1 ). A more detailed analysis, including for instance the im- 
pact of an asymmetric supernova explosion, or the possible 
change in the secondary's mass, is left for further investiga- 
tion. 

The tidal disruption scenario remains a viable, if fine-tuned, 
unified hypot hesis for the S-stars, the runaway and the hyper - 
velocity stars dPeretsll2009HLu et alj2010HZhang et alj2010t) . 
It requires compact binaries on highly elliptical orbits prior 
to their tidal disruption. While binaries forming and migrat- 
ing in a gas disk may naturally become compact (as shown 
in this study), they tend to have small or modest eccentrici- 
ties. Additional eccentricity excitation mechanisms, such as 
resonant relaxation or a hypothetical intermediate-mass black 
hole, need to be invoked. Similarly, the tidal disruption of disk 
binaries gives stars with small inclination relative to the disk 
they formed and migrated in. Again, additional mechanisms 
are required to account for the apparent isotropic distribution 
of the S-stars. In both these respects, the supernova disrup- 
tion scenario requires less restricted boundary conditions. It 
provides a natural way to form bound stars with potentially 
large eccentricity and random inclination (depending on the 
misalignment of the internal and orbital angular momentum 
vectors of the binary prior to its disruption) from circular, 
compact binaries. In the context of the hypervelocity stars, 
the supernova disruption mechanism also provides a natural 
way for a 3 — 4M star to be young at the point of ejection, 
whereas the tidal disruption scenario does not select in favor 
of the young B stars. 
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